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We study the synchronization in a one dimensional array of point Josephson junctions coupled to a common 
capacitor, which establishes a long-range interaction between junctions and synchronizes them. The stability 
diagram of synchronization in a noise-free system is obtained. The current when junctions transform from 
resistive state into zero-voltage state, is then calculated and its dependence on the shunt parameters and the dis- 
sipation of junctions is revealed. In the presence of thermal noise, the synchronized oscillations are destroyed at 
a critical temperature and the system undergoes a continuous phase transition of desynchronization. A possible 
stability diagram of the synchronized oscillations with respect to thermal noise, current, dissipations and shunt 
capacitance is then constructed. Finally we investigate the dynamic relaxation from random oscillations into 
synchronized state. The relaxation time increases with the system size and temperature, but is reduced by the 
shunt capacitor. 
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I. INTRODUCTION 

Josephson junctions are building elements of many elec- 
tronic and electromagnetic devices as well as a candidate for 
quantum computers[jJ. In practical applications, one usually 
integrates large arrays of junctions on a chip to enhance the 
performance, thus coherent operations in these junctions are 
crucial. The synchronization between junctions can be re- 
alized by coupling them to a common resonator, most fre- 
quently through electromagnetic coupling. The common res- 
onator establishes long-range interaction between junctions, 
which then synchronizes them under appropriate condition. 
The junctions arrays have become an extremely important 
playground to understand the synchronization mechanism for 
large population of nonlinear oscillators, partially because of 
the relatively easy experimental realization. 121- fTOl . 

The successful observations of coherent emission from 
cuprate superconductors renew the interests in understand- 
ing the synchronization of arrays of Josephson junctions [llT - 
[TTl . Cuprate superconductors, such as Bi2Sr2CaCu208+5 
(BSCCO), are a natural realization of a stack of Josephson 
junctions of atomic thickness ifTSl \T9 ]. known as intrinsic 
Josephson junctions (IJJs). Because of the large supercoduct- 
ing energy gap, these build-in Josephson junctions can be op- 
erated at frequencies in the terahertz region, where the elec- 
tromagnetic waves have wide applications Il20ll2n . 

Radiation from IJJs occurs in the resistive state. Such a state 
is reached by increasing the bias current above the Josephson 
critical cuiTent and then diminishing it down to the voltage V 
corresponding to the target frequency according to the Joseph- 
son relation w = 2eV/Ji. The resistive state is preserved down 
to the retrapping current below which the system undergoes 
transition into the zero-voltage state. Such a procedure is pos- 
sible because the resistivity of IJJs is very large, i.e. junctions 
are strongly underdamped. Thus the hysteretic behavior al- 
lows us, in principle, to reach a quite low voltage of the order 
of that corresponding to the Josephson frequency (~ 0.1 THz 
for BSCCO). In the resistive state, Josephson plasma of com- 



posite oscillations of Cooper pairs and electromagnetic waves 
is excited. If the plasma oscillations are in-phase, then the to- 
tal radiation power is proportional to the number of junctions 
squared. Below a threshold current called retrapping current, 
the resistive state becomes unstable and the system switches 
into zero-voltage one. Important questions to be addressed 
are: 

1 . What is the retrapping current in the array of point junc- 
tions and how a shunt affects it. 

2. In what parameter region of junction and shunt where 
oscillations of junctions remain synchronized in the re- 
sistive state. 

The stability of synchronized oscillations depends crucially 
on interaction between junctions. The junctions in cuprate su- 
perconductors interact with each other through nearest neigh- 
bor coupling, either inductive or capacitive. These short-range 
interaction however is insufficient to establish a global phase 
coherence Il22ll23l . There are two methods to achieve global 
synchronization by coupling all junctions to a common res- 
onator 

In the first approach, the cavity formed by the superconduc- 
tors single crystal plays a role of the resonator IIT2llT3l . The 
synchronization is realized by the excitation of cavity mode 
in the crystal ll24l[25l . Alternatively, the synchronization can 
be achieved by the radiation fields 1 26 1 and/or by a shunted 
circuit 1 27]. The synchronization by a shunted circuit attracts 
considerable interests, because it can be implemented easily. 

Real junctions involve thermal noise, especially for those 
in high-Tc superconductors. Generally, one expects thermal 
noise broaden the linewidth of oscillating spectrum, or even 
destroys the coherence. It is preferable to have robust coherent 
oscillation against noise. To this end, it is important to know 
how thermal fluctuations destroy the synchronization. 

The dynamical process of building up the synchronization 
is also important for both applications and theoretical un- 
derstandings. For an initial condition that is very close to 
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FIG. 1. (Color online), (a) Schematic view of an array of Josephson 
junctions shunted with a capacitor. The junctions are biased by a dc 
current Ig- (b) The Josephson junction is modeled as a shunt circuit 
of a capacitor, a resistor and a nonlinear Josephson current. 



the fully synchronized state, the relaxation to synchronized 
state can be analyzed based on the standard local stability 
analvsis li26l 1281 lZ9ll . However, for a complete random initial 
state, the dynamic process is highly nontrivial. The system 
may even not relax into the synchronized state. Two ques- 
tions naturally arise: how to reach the synchronized state in a 
controlled way and what is the relaxation time? 

In this paper, we consider a one dimensional array of point 
Josephson junctions coupled to a common circuit. First we 
provide analytical and numerical study on the stability of the 
synchronized state and map out the stability phase diagram. 
Based on the diagram, we derive the dependence of the re- 
trapping current on the shunt circuit. Then we introduce ther- 
mal noise into the system and describe the effect on the syn- 
chronization. Mean-field critical behaviors are identified at 
the desynchronization transition, i.e. transition from the the 
synchronized state to the state with random or partially ran- 
dom oscillations. We reveal the dependence of the transition 
temperature on the shunt circuit and bias cuiTent. Based on 
these results, a possible stability diagram of the synchronized 
oscillations is constructed taking thermal noise into account. 
Finally, we study the relaxation dynamics starting from a dis- 
ordered state, where junctions oscillate randomly. 

The remaining part of the paper is organized as follows. In 
Sec. II, we introduce the model. In Sec. Ill, we perform stabil- 
ity analysis of the synchronized oscillations both numeiically 
and analytically. In Sec. IV, We study the desynchronization 
transition of the coherent state and obtain the corresponding 
transition temperature. In Sec. V, we study the dynamic relax- 
ation from disordered initial state into the synchronized state. 
The paper is concluded by a short summary. 



II. MODEL 

AiTays of Josephson junctions coupled to a common load 
have been extensively studied decades ago 013121, not only 
for their importance for the application in electronic device, 
but also as a fruitful platform to understand the underlying 
synchronization mechanism. These models although are less 
transparent than the well known Kuramoto model 130 1, can 
be realized experimentallv lSTl [32l much easier than the Ku- 
ramoto model. The latter has been realized experimentally 
only very recently |33 1, long after its proposal. Some specific 
configuration of the array, such as one dimensional array of 
Josephson junctions shunted by a serial RLC circuit, can be 
mapped into the Kuramoto model H . 

We consider a stack of IJJs with lateral sizes of order of 
several micrometers. This geometry of junctions is an al- 
ternative route to strong emissions 1 26 1 and has attracted lots 
of attention recently. In this case, the variation of supercon- 
ductivity phase in the lateral direction is small and the junc- 
tion can be approximated as a point junction. The inductive 
coupling |'34'-'3Fl between junctions then vanishes under this 
approximation. Meanwhile, the capacitive coupling ||37l l38l 
between junctions is weak and short-range, thus it can be ne- 
glected in comparison with the long-range interaction medi- 
ated by the shunt circuit. Under these simplifications, a stack 
of IJJs reduce to a serial array of point junctions. 

We study a serial array of point Josephson junctions 
shunted by a lumped C circuit, which is shown schematically 
in Fig. [T] Each junction is modeled as resistively and capac- 
itively shunted circuit. The total current across the junction 
is 
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where V - j^cpk is the voltage of the junction according to 
the ac Josephson relation. Here (pk is the gauge invariant su- 
perconductivity phase difference of fc-th junction, and Rj, Cj 
and Ic are the resistance, capacitance and critical cuiTent of 
the junction respectively. Using the Kirchhoff 's loop law, we 
obtain the equation of motion 



h ^ Q + Icsin(pi^ + 



V = 



2e 



Q 



(2) 



(3) 



where Q is the charge on the shunted capacitance, C, is the 
shunted capacitance and /g is the bias DC current. We have 
introduced the Nyquist noise (white noise) current I^, 

</^> = 0, {q(t)Il,(t')}^(4kBT/Rj)6(t-t')5(k-k'), (4) 

where ks is the Boltzmann constant and T is the temperature. 
We have also assumed that junctions are identical. In the pres- 
ence of the common circuit, small spread in the junction's pa- 
rameters will not destroy the coherent oscillations. 
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We will use dimensionless quantities in the following cal- 
culations. The time is in units of Josephson plasma frequency 
(Op - -\/2eIi./TiCj, current in units of 7^,, capacitance in units 
of Cj, resistance in units of Rj. We then arrive at the dimen- 
sionless version of Eqs. ([3]) and Q 



Ib = Q + sin(;A^ + f3^k +'<Pk+ 



y 



{Il{t)Il,{t')) = ipTSit - t')6(k - k'\ 



(5) 



(6) 



(7) 



where l/yS^ with jS = y/n/(Rj ^/2eI^) is the McCumber 
number which determines the hysteretic behavior of junctions. 
Upon increasing the bias current, the system remains zero- 
voltage until the bias current exceeds the critical current. Then 
the junctions switch into resistive state. The system keeps re- 
sistive even when the bias current is reduced below the crit- 
ical current for a junction with small /3(39]. In the resistive 
state, the superconductivity phase (p^ is rotating accompanied 
by small oscillations. The angular velocity of the rotation for 
ipi; is the same for all junctions determined by the voltage, but 
the phase ipi^ may vary from junction to junction. We will con- 
sider the synchronization of the phase of the junction arrays 
in the following. 

In this model, all junctions are coupled to the capaci- 
tor C.5, which establishes mutual interaction among all junc- 
tions. Thus the effective dimensionality of the system is infi- 
nite, which is crucial for the synchronization ll22l |23]| . An- 
other consequence of this mean-field behavior is the per- 
mutation symmetry, i.e., all junctions are biased by the 
same external current and the current in the shunt circuit. 
The exchange of any pair of junctions in the circuit does 
not change the topology of the circuit. If configuration 
((^1, 02, 0/, is a solution to Eqs. (|5]l, (j6|l and 
(jTji, then ((f)i,<p2, ■■■,</>!, ■■■,<pi, ■■■,<Pn) is also a solution. This 
symmetry greatly simplifies the stability analysis as will be 
shown below. 

Apparently, Eqs. (|5]l, (|6]l and (|7]i always have a trivial so- 
lution that all junctions oscillate out of phase, and the dynam- 
ics of each junction is independent because the current in the 
shunt circuit vanishes. However suppose at some instance, 
a small population of junctions oscillate with the same phase, 
then the capacitance C, acquires energy, which is proportional 
to the number of in-phase junctions squared. Now the capac- 
itance is able to attract more junctions to oscillate at its phase 
and in turn its energy increases further. This is a positive feed- 
back process with explosive increases of energy in the capac- 
itance and an avalanche of junctions oscillating coherently. 
Therefore we expect in certain parameter space, the out-of- 
phase oscillations lose stability and synchronization sets in. 
The qualitative picture will be elaborated in subsequent sec- 
tions. 



III. STABILITY OF THE SYNCHRONIZED STATE 

In this section, we analyze the local stability of the coherent 
oscillations of all junctions in the absence of thermal fluctua- 
tions, T -Q. The local stability is determined by the dynamics 
of the system in the vicinity of the trajectory of the uniform 
solution. We consider the uniform oscillations (pk = (f>o, where 



(A^C, + 1)00 +y60o + sin 00 = Ig, 



(8) 



with being the number of junctions. The junction coupling 
strength is enhanced by a factor of A^, in accordance with 
the typical behavior in the mean-field theory. We then add 
small perturbations 6^ to the uniform solution and determine 
the time evolution of the perturbations. The equations for the 
perturbations read 



'Sk + 136k + cos(0o)5a. + C, ^ 6i = 0. 



(9) 



(=1 



The permutation symmetry between junctions allows us to de- 
couple Eqs. (|9]l by introducing the quantities = 6k+\ - dk 

and cr - h Yi 6k- We obtain equations for = A and cr: 



A + ySA + cos(0o)A = 0, 



& + P& + cos(0o)cr + CsNa- - 0. 



(10) 



(11) 



If A diverges with time, the uniform solution becomes unsta- 
ble. On the other hand, if cr diverges while A decays with 
time, the synchronization is kept and the system transits into 
another synchronized state if it exists. We are interested in the 



coherent oscillation, and we will only focus on Eq. (10 1 in 



the later analysis. We will solve Eq. (lOi for weak oscilla 



tions both analytically and numerically based on the Floquet 
theorem. 



A. Analytical treatment 

We consider the region where the amplitude of Josephson 
oscillation is small. The solution of Eq. (|8]l in linear approxi- 
mation can be written as 



ia-tot + A expiicot) 



with 



« 1. 



(12) 



(13) 



The frequency w is determined by the DC current conserva- 
tion 



/b =ySw-HRe[A]/2. 



(14) 



Substituting Eq. ( 12 1 into Eq. ( 10 1, we get the equation for A 



A+/3A + 



^[e''^' + e-''^')- j.[e^''^' -1)a 



A = 0. (15) 
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FIG. 2. (Color online). The largest Floquet exponent calculated by 
detD = with D given by Eq. ijTSj [red line], and by numerical 
calculation using the Floquet theory in Eq. (21) [symbols]. For the 
exponent smaller than 0, the uniform oscillations are stable. 



The coupling of perturbations to the oscillation exp(iojt) in- 
duces higher frequency harmonics. The general solution for 
A is 



FIG. 3. (Color online). Stability diagram of the uniform solution 
in the absence thermal fluctuations. Light blue/pink/orange region 
denotes complete synchronization/zero-voltage state/partial synchro- 
nization. The blue line is the stability boundary of the uniform solu- 
tion calculated by the Floquet theory and the open red circle is the 
retrapping current determined by direct calculations of Eqs. l |5|7^ . 
The dashed line the retrapping current determined by Eq. l[8j while 
the dotted line is the retrapping current for a single junction. Here 
C, = 3/N. 



k=-oo 



Jk(i)t 



(16) 



The stability is determined by the spectrum of perturbations 
Q. In the framework of the Floquet theory 1401 . we call Im(Q) 
as the Floquet exponent. The uniform solution is stable if and 
only if the largest Floquet exponent is negative, i.e. Im(Q) < 
0. One may easily identify Im(Q) as a relaxation time and 
Re(Q) as an energy gap of perturbations. 

To obtain Q, we plug Eq. ( 16 1 into Eq. ( 15 i and compare 



each frequency component. Then we have the following linear 
equations for the coefficients a^. 

1 A 
-(kLo-D.fak+i(k(Aj-Q.)/3ak+- (flu + flt+i)-— (fl^-a - flt) = 

2 2i 

(17) 

The existence of nonzero solution of requires the determi- 
nant of the coefficient matrix vanishes, det D = with 
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(18) 



with solutions 



Q 



2C,N 1 
(C,iV+l) w2 



(20) 



We see that the uniform oscillations are always stable for non- 
zero Cj in the region of w » 1. In the limit a> ~ Ib//3 — > oo, 
the largest Floquet exponent approaches zero, and the solution 
becomes neutrally stable. 

Near the stability boundary where the largest Im(Q) 
changes sign, one has to keep higher harmonics in Eq 
because cj ~ I 
expansion in Eq 



(16l 



But for CsN » 1, one can still use the linear 
Under these conditions, the stability 



12i 



boundary can be determined by the numerical calculation of 
detD = 0. 



B. The Floquet theory 



Equation ( 10 1 can also be interpreted as a particle moving in 
a periodic potential with period T. Then we can apply the Flo- 
quet theorem (Bloch theorem) to extract the exponents. The 
solution has the form 



andQ = -(kcj-Q.)^ + i(ka>-Q.)/3 + A/{2i). The solution gives 
the spectrum of the perturbation. 

In the region of w » 1, the frequency modes with k -Q,+l 
are dominant and higher harmonics may be truncated. We 
obtain a second order equation for Q 



+ ipO. 



A 1 



C,N 



1 



(C,N + 1) 2w2 ' 



(19) 



A(f) = exp(^if)3;i(f) + exp(A2t)y2(t), 



(21) 



where yi (t) and y2(t) are periodic functions with period T, and 
the exponents Ai and A2 follow Ai + A2 - -p according to the 
Floquet theorem. il40l When no dissipation is present yS = 0, 
the dynamics is time reversal and Ai + A2 = 0. When the 
dissipation is involved, the volume of phase space is shrinking 
with a rate /3, thus the two exponents follow Ai + A2 = -/?. 
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The exponents can be computed numerically as follows. 
We first calculate the trajectory of 00 in Eq. ([8]). Then we 
calculate two trajectories of Aa(f) and Afc(f) with two differ- 
ent initial conditions Afl(fo) - 0, Afl(fo) = 1 and A/,(fo) - 1, 
Ai(fo) = 0. These two trajectories obey[:40l 



AJt + T),Ah(t + T) 
Aa(t + T),Ah{t + T) 



)-<^'(i:!;!:i:;;!)' 



with F being a coefficient matrix, which can be evaluated by 
inverting Eq. (22i because the trajectories of A^, and A;, are 
known. Ai and A2 are just the eigenvalues of the matrix F. 



We have compared the results obtained by Eq. (22 1 and those 
by analytical calculations. Both methods give the consistent 
results as shown in Fig. |2] 
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C. Stability diagram 

The stability analysis above does not tell us what the fi- 
nal state is when the uniform solution becomes unstable. To 
answer this question, we solve Eq. (jSjl and (j6jl with I'^ = 
directly by numerical simulation. The stability diagram then 
is constructed, and is depicted in Fig. [3] For a sufficient large 
Ib thus <y » 1, the uniform solution is stable as described by 



Eq. ( 20 1. For a small /3, upon decreasing Ig, the uniform oscil- 
lations become unstable below the retrapping current /,. < Ic 
and the system evolves into zero-voltage state. For a large 
yS, the uniform oscillation loses stability at /, > Ic, where no 
zero-voltage state is available for the system to go. In this 
case, the system becomes partially synchronized with a frac- 
tion of junctions oscillating in-phase, while the others do out- 
of-phase oscillation. When Ib is reduced further, the partial 
synchronization becomes unstable and the system is retrapped 
into zero-voltage state at Ib - Ic- 

Let us discuss the transition from the complete synchro- 
nization to the partial synchronization when /3 is large. To 
characterize the partially synchronized state, we introduce the 




FIG. 5. (Color online). Dependence of the boundary retrapping cur- 
rent Ir on the shunt capacitance for several p's. Lines are obtained 
with the Floquet theory and symbols are direction simulations of Eqs. 
l |5|7[ ). The region above the line corresponds to the uniform oscilla- 
tions while zero-voltage state below the line. Here A' = 200. 



order parameter which is widely used in literatures ll30l 



1 ^ 



rit)exp[ie(t)] = - 2^exp(i(Pj) 



(23) 



Here r is positively defined. We compute the average of r(f) 



(r) 



1 r'l 
Jo 



dt r(t). 



(24) 



and take fy +oo. 

The IV and the corresponding order parameter are shown in 
Fig. |4] When the complete synchronization becomes unsta- 
ble, a sharp jump of voltage is observed, associated with de- 
crease of the order parameter. The reduction of voltage when 
the system becomes partially synchronized can be understood 
as follows. At a given voltage, the shunt capacitor reduces 
the plasma oscillation amplitude depending on the number of 
the synchronized junctions, as described by Eq. ( [T3] l. For 
the uniform oscillations, the suppression is largest and the DC 
current induced by the Josephson oscillation is reduced sig- 



nificantly according to Eq. ( 14 1. For the partial synchronized 



oscillations, the DC current is larger than that of the uniform 
oscillations. Therefore when one biases the array with a fixed 
current, the voltage of the uniform state increases compared 
with that in partial synchronized state. 
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FIG. 4. (Color online). IV curve(black) and dependence of the order 
parameter on the voltage(red) for p = 2.0 and Cj = 3/A'. 



D. Retrapping current 



According to Eqs. (13 i, ( 14 1 the amplitude of Josephson os- 
cillation increases with decreasing Ib, To achieve the strongest 
oscillation, one would like to know how small current one can 
achieve in order to support the resistive state. 
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FIG. 6. (Color online). IV curves for (a) jS = 0.02 and (b) /3 = 1.0 for 
several typical values of C,. Here A' = 200. 



For a single junction, the dynamics is equivalent to a parti- 
cle sliding down in the damped inclined washboard potential. 
It shows hysteretic behavior for small /3, i.e. the system re- 
mains resistive even Is < h- The system evolves into the su- 
perconduction state at a current 1^ > 0, where the input power 
is insufficient for the phase particle to move in the damped 
tilted washboard potential. The retrapping current for a weak 
damping is given by 1^ ~ 1.48y6.|4r| On the other side, the 
dynamics becomes overdamped for a large /3, and the system 
comes back to zero-voltage state once Is < h- The depen- 
dence of Ir on y6 for a single junction is shown by dotted curve 
in Fig. [3] 

For junction array shown in Fig. [T] if the uniform solution 
is always stable in the whole cuiTent region, the retrapping 
cuiTent will be the same as in a single junction case with an 
effective - /?/ '\/CsN + 1 normalized by the shunt capacitor 
(dashed line in Fig. |3]l. In fact, the uniform solution loses 
stability at Is (blue line in Fig. [3]l and the system evolves into 
the zero-voltage state. Therefore Is is the genuine retrapping 
cuiTent for the present junctions aiTay, and can be measured 
experimentally. 

How to decrease or is it possible to shift the stability 
boundary in Fig. [Sjleftward? One recalls that the shunt capac- 
itor induces interaction between junctions. By increasing the 
coupling constant Cs, one would expect that the stable region 
enlarges and the stability boundary shift leftward. We study 
the dependence of Is on Cs, and the results are presented in 
Fig. |5] ForyS > 0.5, the retrapping cuiTent decreases with Cs, 
while it increases with Cs for smaller /3. A qualitative picture 
for this unexpected non-mo no tonic dependence is as follows. 
Equations 



and (10 1 with - can also describe the 
stability of the resistive state for a single junction, where the 
retrapping current is enhanced by the dissipation as shown 
by the dotted line in Fig. |3] We rewrite Eq. (|8]l to a 
form equivalent to the single junction case by rescaling the 
time t «- y/NCs 
PI iNCs + 1 
then acquire a form 



+ It', and with an reduced dissipation /3' 
The dynamics of the perturbations Eq. (10 1 



1 



NC, + 1 



A + j8'A + cos(0o)A - 0. 



(25) 
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FIG. 7. (Color online). Same as Fig. |5]but with a shunt RC circuit. 
Parameters are indicated in the figures. 



of the mass makes the system more vulnerable to the oscil- 
latory potential cos(0o)A^/2, therefore tends to increase the 
retrapping current. For the junction array with weak damping, 
the effect of the reduction of mass dominates and the retrap- 
ping current is increased by the shunt capacitor. On the other 
hand, for the junction array with strong damping, the reduc- 
tion of the damping by the shunt capacitor dominates and the 
retrapping current is decreased. 

The IV curves for several typical values of C, are shown in 
Fig. |6] As seen in Fig. |6], the IV curves with C, = de- 
viate from the asymptotic linear behavior Ig ~ w/yS at small 
o), which indicates strong plasma oscillations according to Eq. 
( 14 1. For j6 = 0.02 the IV curves in Fig. [6ja) behave differ- 
ently near the trapping point for C, - and for nonzero Cs, 
where in the latter case the IV is linear down to the retrapping 
current. This linear dependence near the retrapping point with 
(!) ~ 1 is due to the suppression of the oscillation amplitude 
by the shunt capacitor Cs- However, for y6 = 1.0 the IV curves 
remain nonlinear near the retrapping point even for the same 
Cs- For a large y6, the retrapping voltage is much smaller than 
that of small yS, thus results in stronger plasma oscillations. 



Thus the presence of shunt capacitor reduces both the effective 
mass of perturbations and damping coefficient. The reduction 



E. RC circuit 

In real devices, the shunt circuit also carries finite resis- 
tance. We perform similar analysis of the stability diagram 
and the retrapping current, when a resistor with resistance R 
is serially connected to the capacitor in the load circuit. A new 
term QR should be added to the right-hand side of Eq. (j6]l. 

The resulting stability diagram is qualitative the same as 
Fig. [3] For a given R, the retrapping current increases with Cs 
for a small /3 while decreases for large /3, as depicted in Fig. [7] 
Comparing Fig. [Tjwith Fig. [5] we can see that the shunt resis- 
tor increases the retrapping current, because the dissipation of 
the system is increased by the resistor. 
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FIG. 8. (Color online). Stability of the uniform oscillations when an 
array of Josephson junctions are shunted by an LRC circuit. The 
region below lines are stable. 



FIG. 9. (Color online). Dependence of the order parameter on tem- 
perature with different system sizes. Inset is a double-logarithm plot 
of the reduce temperature T,„ - T and order parameter. Here /g = 1.5, 
/3= 1.0 and C, = 3.0/N. 



F. LRC circuit 

Another interesting case is an array of Josephson junctions 
shunted by a LRC circuit, which introduces a characteristic 
frequency oji. - 1/ VAsCs- The LRC circuit can represent the 
cavity intrinsically formed by the single crystal of BSCCO 

ma [m nam. 

The stability of the uniform oscillations can be obtained 
similarly. We consider the case with w » 1 and y6 «: 1 so 
the analysis in Sec. III(A) is applicable. The dynamics for 



small perturbations is still given by Eqs. ( 10 1 and ( 12 1 with a 
modified amplitude 



(26) 



Nop- 



L,(a)J-a)2)+(ji7f 



The stability is determined by Eq. (19i and the results are 
shown in Fig. |8] When co <«c Uc, the LRC circuit behaves 
as a RC circuit, which always stabilizes the uniform solution 
for a small R as given by Eq. (20i. On the other hand, for 



Lo » LOc, the LRC circuit behaves as a LR circuit which makes 
the uniform solution unstable for a large Lj. When u) ~ a>c, 
the stability depends on the quality factor R of the LRC cir- 
cuit. Small quality factor (large R) makes the synchronization 
difficult. 



IV. EFFECT OF THERMAL NOISE 

Real circuits inevitably involve noise because of resistiv- 
ity caused by qu as ip articles. This leads to diffusive dynamics 
in the phase space and destroys the synchronization at a cer- 
tain critical point. To study the effect of noise, knowledge 
of attractors in the phase space is necessary. An attractor at- 
tracts trajectories nearby and the volume of phase space that 



the attractor attracts defines the basin of attraction of the at- 
tractor The phase space is covered by basins of attraction 
and the boundary between basin of attraction is called separa- 
trix. For a small /3, we have already identified two attractors 
with one being the zero-voltage state, and the other uniform- 
oscillation state. For a large /3, additional attractor with partial 
synchronous oscillations appears. 

Suppose the system initially at the uniform-oscillation state, 
and then we turn on thermal noise. The noise perturb the sys- 
tem away from the uniform state. However, the deviation from 
the attractor is penalized by the action S. The comparison 
of the action S for the system moving from the attractor to 
the separatrix with the noise strength defines three distinct re- 
gions. 



0.015 - 



0.010 



-N=3200 
-N=1600 




0.05 



FIG. 10. (Color online). Dependence of the fluctuations a-^ on the 
temperature. Here /g = 1.5,/? = 1.0 and Cj = 3.0/N. 
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FIG. 11. (Color online), (al), (bl) and (cl): dependence of T„, on the bias current, shunt capacitance and fi respectively. (a2), (b2) and (c2) 
are the conesponding largest Floquet exponent. Other parameters used are shown in the figure. 



1. kbT <«; 5, in this case the possibility of thermal escape 
is extremely small and this region is described by the 
reaction-rate theory 1 44 1. Especially, when the thermal 
activation between two attractors a and b is asymmetric, 
that is the action S a '» Sh, the system will spend much 
longer time in the attractor b. This is the situation of 
retrapping from resistive state to zero-voltage state for 
small y6 <sc 1 discussed before. When the bias current 
Ib is close to the retrapping cuiTent Ib - Ir ^ 1, the 
resistive state is about to lose stability. So the presence 
of weak noise will destabilize the resistive state and the 
system evolves into the zero-voltage one. On the other 
hand, the energy barrier for the system to transform 
from zero-voltage state to resistive one again is large 
when Ib Ic, and the noise are not strong enough to 
promote such a transition. So the system remains zero- 
voltage. Thus the thermal noise increase the retrapping 
currentll45l. 

2. k[,T » 5, in this region the thermal energy is large 
enough to kick the system off the attractor of the co- 
herent oscillation, and the synchronization is destroyed. 
The temperature at which the synchronization is de- 
stroyed is the synchronization-desynchronization tran- 
sition temperature r„,. 

3. for T < T,„, the uniform oscillations survive. The noise 
current excites perturbations and the system frequently 
deviates from the attractor. The dynamics of the pertur- 



bations are described by Eq. (j9]l. These perturbations 
broaden the linewidth of the frequency spectrum. The 
linewidth at w » 1 can be estimated as follows. For 
w » 1, the IV is linear, so the noise current /" induces 
a noise voltage /"//?■ From the ac Josephson relation 
dt4> = 2eV/ti, one easily obtain that the linewidth in- 
creases linearly with T for Gaussian white noise. 

In the presence of noise, the equations of motion become 
stochastic, and it is natural to describe the dynamics in term 
of a probability density in the phase space. The flow of the 
probability density is governed by the Fokker-Planck equa- 
tion. However analytical calculations of the coupled non- 
linear partial differential Fokker-Planck equation is difficult. 
In this section we will use numerical simulations as a main 
workhorse, and we will also provide qualitative analysis to 
understand the numerical results. We first consider the desyn- 
chronization transition of the synchronous state and the criti- 
cal behavior at the transition. We then find a correlation be- 
tween the largest Floquet exponent and the transition temper- 
ature. Finally a stability diagram of the uniform oscillations 
with respect to noise is constructed. 



A. Synchronization-desynchronization transition 



To study the synchronization-desynchronization transition. 



we evaluate the order parameter defined in Eqs. ( 23 1 and ( 24 1 
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and its standard deviation 



(27) 



which is similar to the susceptibiHty defined in spin systems. 

We solve numerically Eqs. (j5]l, (j6]l and (j7|i, and derive (r) 
and cTr at diff'erent T . The results are presented in Figs. |9] 
and[TO] We also check the finite-size effect with different A^'s. 
The finite size eff'ect is prominent around T,„. The synchro- 
nized oscillations is continuously suppressed by the thermal 
fluctuations. At T,„ the synchronized oscillations become un- 
stable, and the system undergoes a continuous transition into 
random oscillations. Since cr^ serves as a measure of the fluc- 



tuation effect, it reaches maximum at r,„, as shown in Fig. 10 



Practically cTr provides a convenient way to determine r„, es- 
pecially for a small system where the transition is obscured by 
the finite-size eff'ect. 

Once we identify the desynchronization transition as a crit- 
ical phenomenon, we can define the exponent 



(28) 



In Fig. |9] we obtain fSc ~ 1 /2 which is consistent with the 
mean-field theory. 

We then study which factors determine T,„ and how to en- 
hance T,n- As discussed at the beginning of this section, T,„ is 
given by the action for the system moving out of the attractor 
Thus knowledge about the whole basin of attraction is needed. 
However, it is still conceivable that the local slope near the 
attractor may to certain extent reflect the global structure of 
basin. The local slope is just the largest Floquet exponent ob- 
tained in the previous section. Therefore one expects that the 
smaller the exponent, the higher We find numerically that 



it is indeed the case, as shown in Fig. 1 1 



The correlation between the largest Floquet exponent and 
T,„ can be understood in terms of the local stability analysis. 
The small perturbations to the uniform oscillation decay 'q ~ 
exp(/l20 with A2 <Q being the largest Floquet exponent. This 
is equivalent to the relaxation of a particle in the parabolic 
potential d,^ = -dV/dq with V(q) - -/1 2^/2. The slope 
-A2 > measures the depth of the potential. Thus it is more 
robust against noise for a larger -/I2. 




FIG. 12. (Color online). Possible stability diagram of the uniform 
solution (a) at a given C.,, and (b) at a given p. The region inside the 
green surface corresponds to stable uniform oscillations. 



B. Stability pliase diagram with noise 

Based on the previous analysis, we discuss the stability 
phase diagram of the synchronization in the presence of ther- 
mal noise. For a given C.,, when the bias current is increased, 
the system approaches the synchronous state, where the asso- 
ciated Floquet exponent changes from positive to negative at 
the stability boundary. If the current increases further, it reach 
the maximal value -/3/2. For a sufficient large current, the sys- 



tem becomes neutral stable according to Eq. (20i. Therefore 
the critical temperature first increases and then decreases with 
the current. The corresponding stability diagram is shown in 
Fig. 12 a). Meanwhile, the shunt capacitance plays a role 



of coupling strength, so r,„ increases with C,. Keep in mind 
that the current at the stability boundary is the retrapping cur- 
rent. At a given T, the retrapping current increase with C, for 
a small /3 , while it decreases for a large /?. Based on these 
observations, we construct the phase diagram of the coherent 



oscillations for a given /3, which is sketched in Fig. 12 b). The 



region enclosed with green surface represents a stable syn- 
chronization. 

To enhance T,„ a larger /? is helpful since the maximal Flo- 
quet exponent is -y6/2. One should also adjust the current 
accordingly to ensure that the maximum is reached. For a 
given operating frequency, one may increase C, to enhance 
the thermal stability, at sacrifice of the oscillating amplitude. 



V. DYNAMIC RELAXATION 

So far we have concentrated on the stability of the uni- 
form oscillation, and investigated the dynamics of perturba- 
tions around the uniform state. However, in most applications, 
the initial condition cannot be guaranteed as the state of uni- 
form oscillations. For instance, when we ramp up the current 
and bias all junctions in the resistive state, the initial state may 
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T=0 T=0,0005 T=0.002 




FIG. 13. (Color online). Time evolution of the distribution of phase difference Eq. \29\ , starting from completely random state at several 
temperatures. Here fj = 0.02, Ig = 1.5 and = 3/A'. 



be far away from the uniform state in the phase space. There- 
fore it is important to understand how the system approaches 
the uniform state. In the present study, we focus on the relax- 
ation from disordered state (all junctions oscillate out-phase) 
into ordered state (all junctions oscillate uniformly). For a 
system whose final ordered state is in equilibrium, this is a 
phase ordering phenomenon. The kinetics of phase ordering 
has been extensively studied decades ago in spin systems, and 
they can be described by universal scaling behavior|46 1. How- 
ever the relaxation dynamics is not very clear when the final 



state is out of equilibrium. 

To reach the uniform state, the initial state must be in the 
basin of attraction of the uniform state. This can be realized 
by operating all junctions at resistive state. We prepare the ini- 
tial state with arbitrary nonzero (r) <s 1 . We also give initial 
velocity to all junctions, as such the system falls into the basin 
of attraction of the uniform state. Let us first consider dynami- 
cal relaxation obtained by computer simulation. We introduce 
the distribution of the phase difference between junctions 



(29) 



(a) 


















■ T=0 




• T=0.001 




A T=0.002 




1000 0.000 



0.001 0.002 

T 



0.03 0.04 



FIG. 14. (Color online), (a) Dependence of the relaxation time on the 
system size with different temperatures. Symbols are numerics and 
lines are the best fitting, (b) Dependence of the relaxation time on 
the temperature with the system size of N = 200 and with different 
shunt capacitance, (c) Speedup of the relaxation by increasing the 
shunt capacitance. Here the system size is N = 200 and temperature 
T = 0.002. All these results are obtained with /3 = 0.02 and 1b = 
0.15. 



with Aij - 4>i - <pj- The time evolution of P(<p) is depicted in 
Fig. 13 Initially the distribution is flat indicating a disordered 
phase. This flat distribution does not change with time too 
much at the beginning, but then it suddenly becomes sharp. 
Finally it reaches a steady distribution with finite width de- 
pending on the temperature. 

Qualitative picture of the relaxation can be obtained based 
on the local stability analysis presented in Section III. Suppose 
we have a small synchronized cluster of junctions with pop- 
ulation A^, and the rest of junctions oscillate randomly. This 
small cluster serves as a seed of the nucleation, and deliver 
energy into the shunt capacitor, which in turn attracts nearby 
out-phase oscillators into the cluster. The growth rate of the 
synchronized population can be estimated by the local stabil- 
ity analysis by replacing with «(f), where n(f) is size of the 
cluster at time t. The time evolution of the population of the 
cluster follows 



n(t + dt) - n{t)exp{-Adt) « «(f)(l - A{n)dt), 



(30) 



where dt is a small time step and A(n) < is the largest FIo- 
quet exponent with cluster size of n. Then the time required 
for the system achieves global synchronization is 



dn . 

/V nA{n) 



(31) 



Several observations are in order First, the synchronization 
time increases with the total number of junctions A^. Sec- 
ondly, since A{n) < decreases monotonically with n and 



11 



then saturate at -yS/2 [see Eq. (20 1], the initial relaxation is 
slow and it gradually speeds up, in accordance with Fig. [T3] 
Thirdly, in the presence of thermal fluctuations, thermal noise 
may kick oscillators out of the synchronized cluster Thus the 
increase rate is reduced and relaxation time increases. 

To quantify the relaxation process, we define the linear re- 
laxation functionf??! 



A(f) = 



<r(0) - <r(cx,)) 
<r(0))-<r(cx,))- 



(32) 



It starts from unity at f = and decays to in the steady state. 
The relaxation time is defined as 



oo 

/ 



A(t)dt. 



(33) 



For an exponential decay, the definition above is equivalent 
to the conventionally defined relaxation time. Two-stage re- 
laxation for A(t) is found at T = 0. First r increase from 
to a value close to 1, where the local stability theory applies. 
Then the system relaxes into the ordered state exponentially 
with the exponent given by the Floquet exponents. Thermal 
fluctuations smear the distinction of the two-stage relaxation. 

We numerically calculate A(t) and compute r. The depen- 
dence of T on the number of junctions A^, temperature and 
shunt capacitor is plotted in Fig. [14] t increases with as ex- 
pected from the qualitative estimate above. Furthermore the 
relaxation time follows a power law r ~ A^"' . The exponent 
z'iT) increases with T. t increases with temperature. At T„, 
it diverges and then drops, (the relaxation time above T,„ is 
not very meaningful because the final state is also disordered.) 
Critical behavior is also identified for t near r„,. On the other 
hand, t decreases with Cj which suggests a practical way to 
speed up the relaxation. This can be explained by regarding 
Cj as a coupling strength of the system. A larger therefore 
increases the rigidity of the uniform solution. 



VI. CONCLUSION 

In short, we have studied the synchronization of one di- 
mensional array of point Josephson junctions coupled to a 



shunt capacitor In the case of noise-free system, a stability 
phase diagram of the uniform oscillation is constructed. For 
strong damping, after the uniform solution becomes unstable, 
the system evolves into partially synchronized state. When the 
bias current is reduced below the Josephson critical current, 
the system becomes zero-voltage. For weak damping or mod- 
erate damping, after the instability of the uniform solution, the 
system evolves into the zero-voltage state. At transition the 
current is the experimentally measurable retrapping current. 
The retrapping current is increased by the shunt capacitor for 
weak damping (/3 < 0.5), while it decreases for moderate and 
strong damping (fS > 0.5). Thus transport measurement pro- 
vides a convenient probe of the underlying dynamics. Similar 
results are obtained when a resistor is serially connected to the 
shunt capacitor. 

In the presence of strong thermal noise, the coherent os- 
cillation is destroyed through a second order phase transition. 
The critical exponent for the order parameter is 1 /2 in accor- 
dance with the mean-field theory. We also find the fluctuations 
of the order parameter r showing a maximum at the transition, 
which may serve as a convenient quantity to locate the transi- 
tion temperature. For a smaller relaxation time in the case of 
weak perturbations, the transition temperature is higher The 
results suggest several possible ways to enhance the thermal 
stability. 

The dynamic relaxation from a disordered phase to ordered 
state is then investigated. The relaxation time increases with 
the system size by a power law. It also increases when the sys- 
tem approaches the transition temperature from below. One 
may speed up the relaxation with a larger shunt capacitance. 

Finally, a possible phase diagram of the uniform solution is 
proposed when thermal fluctuations are involved. Our results 
are of importance for the design of useful superconducting 
devices based on Josephson junctions arrays. 
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